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Finding the first time a fluctuating quantity reaches a given boundary 
is a deceptively simple-looking problem of vast practical importance 
in physics, biology, chemistry, neuroscience, economics and industry. 
Problems in which the bound to be traversed is itself a fluctuating 
function of time include widely studied settings in neural coding, 
such as neuronal integrators with irregular inputs and internal noise. 
We show that the probability p(t) that a Gauss- Markov process will 
first exceed the boundary at time t suffers a phase transition as 
a function of the roughness of the boundary, as measured by its 
Holder exponent H, with critical value H c = 1/2. For smoother 
boundaries, H > 1/2, the probability density is a continuous func- 
tion of time. For rougher boundaries, H < 1/2, the probability is 
concentrated on a Cantor-like set of zero measure: the probability 
density becomes divergent, almost everywhere either zero or infin- 
ity. The critical point H c = 1/2 corresponds to a widely-studied case 
in the theory of neural coding, where the external input integrated 
by a model neuron is a white-noise process, such as uncorrelated 
but precisely balanced excitatory and inhibitory inputs. We argue 
this transition corresponds to a sharp boundary between rate codes, 
in which the neural firing probability varies smoothly, and temporal 
codes, in which the neuron fires at sharply-defined times regardless 
of the intensity of internal noise. 

random walk | first-passage time | phase transition | neural code 

Abbreviations: FPT, first passage time; OUP, Ornstein-Uhlenbeck process; LIF, leaky 
integrate and fire neuron 

A Brownian process W(t) which starts at t = from 
W{t = 0) =0 will fluctuate up and down, eventually 
crossing the value 1 infinitely many times: for any given real- 
ization of the process W there will be infinitely many different 
values of t for which W(t) — 1. Finding the very first such 
time, 

r = inf{£| W(t) = 1} 

known as the first passage of the process through the boundary 
L = 1, is easier said than done, one of those classical problems 
whose concise statements conceal their difficulty [TJ [21 EJ [4] . 
For general fluctuating random processes the first passage 
time problem (FPTP) is both extremely difficult [5j EJ [7J El E] 
and highly relevant, due to its manifold practical applications: 
it models phenomena as diverse as the onset of chemical reac- 
tions 10, 11, 12, 13, 14 , transitions of macromolecular assem- 
blies [I!l[Il[T3mil9], time-to- failure of a device [20] EQ E2] , 
accumulation of evidence in neural decision-making circuits 
[23] . the "gambler's ruin" problem in game theory [24], species 
extinction probabilities in ecology [25], survival probabilities 
of patients and disease progression 26 , 27, 28 , triggering of 
orders in the stock market [29, 30, 31 , and firing of neural 
action potentials [32 , 33, 34, 35 , 36, 37]. 

Much attention has been devoted to two extensions of this 
basic problem. One is the first passage through a stationary 
boundary within a complex spatial geometry, such as diffu- 
sion in porous media or complex networks, as this describes 
foraging search patterns in ecology 38 , 39 , and the speed at 



which a node can receive and relax information in a complex 
network [40] [41] . 

The second extension is the first passage through a bound- 
ary that is a fluctuating function of time 42 , 43] [44], a problem 
with direct application to the modeling of neural encoding of 
information 45 , 46 . This problem and its application are the 
subject of this paper. The connection arises as follows. The 
membrane voltage of a neuron fluctuates in response both 
to synaptic inputs as well as internal noise. As soon as a 
threshold voltage is exceeded, nonlinear avalanche processes 
are awakened which cause the neuron to generate an action 
potential or spike. Therefore the generation of an action po- 
tential by a neuron involves the first passage of the fluctuating 
membrane voltage through the threshold. This dynamics of 
spike generation underlies neural coding: neurons communi- 
cate information through their electrical spiking, and the func- 
tional relation between the information being encoded and the 
spikes is called a neural code. Two important classes of neural 
code are the rate codes, in which information is only encoded 
in the average number of spikes per unit of time (rate) with- 
out regard to their precise temporal pattern, and the temporal 
codes, in which the precise timing of action potentials, either 
absolute or relative to one another, conveys information. 

Central to the distinction between rate and temporal codes 
is the notion of jitter or temporal reliability. This notion orig- 
inates from repeating an input again and again and aligning 
the resulting spikes to the onset of the stimulus. Time jitter- 
ing is assessed graphically through a raster plot and quanti- 
tatively in a temporal histogram (PSTH) which permits veri- 
fying the temporal accuracy with which the neuronal process 
repeats action potentials. A fundamental observation is that 
the very same neuron may lock onto fast features of a stimulus 
yet show great variability when presented with a featureless, 
smooth stimulus [33] . These two are extreme examples from a 
continuum — the jitter in spike times depends directly on the 
stimulus being presented [47] . 



First passage through a rough boundary 

We shall make use of a simple geometrical construction, map- 
ping the dynamics of a neuron with an input, internal noise 
and a constant threshold voltage, onto a neuron with inter- 
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nal noise and a fluctuating threshold voltage; the construction 
thus maps the input onto fluctuations of the threshold. We 
use as our model neuron the leaky integrate- and- fire neuron 
(LIF) , a simple yet widely-used [36l l47l ft51 l49l f50l f5TT f52l f53l 
E2 El [561 [53 Ell [59] model of neuronal function defined by 



V 



-aV + I(t) + 



[1] 



where V is the membrane voltage, 1/a is a decay time given 
by the RC constant of the membrane, I the current that the 
neuron receives as an input through synapses, and £ an inter- 
nal noise. When V first reaches a threshold value I an action 
potential is generated, and the voltage is reset to zero. The 
nonlinearity of the model is concentrated on the spike gen- 
eration and subsequent reset, so that between spikes we can 
integrate separately the effect of the input and of the noise: 

V = Vi + Vz Vi = -aVi + I(t) V$ = -aVz + £t 

Because the input I(t) is fixed, the Vi equation needs to be 
solved just once. Then the problem of V(t) reaching the 
threshold I can be recast as V% reaching the boundary I — Vi: 
we have transformed a problem with an input and a con- 
stant threshold into a problem with no input and a fluctuating 
threshold l — Vi. The reset operation V = I — >> V = becomes 



V* 



V T 



-Vi (see Appendix) 



These considerations lead us to examine the problem of 
the first passage time through a fluctuating threshold. In or- 
der to develop some intuition about the problem, we are going 
to break it up into two parts, a "geometrical optics" part, in 
which most first passages can be accounted for by simple "vis- 
ibility" considerations, and a "diffractive" correction in which 
we take into account that random walkers can turn around 
corners. The geometrical part is simple: most first passages 




Fig. 1. How a random walk V first hits a moving boundary L. In all panels, 
time t is horizontal, the process V and the boundary L vertical. (A) It is highly 
probable to hit the left flank of a minimum, as the walkers are moving left to right 
and from the bottom up. (B) Each minimum "casts a shadow" behind it, so that 
hitting some features behind may be hard, as it requires missing the minimum, then 
rising sufficiently high to hit the second feature. (C) Hitting the right (rising) flank of 
a minimum is hardest, since it requires missing the minimum narrowly, then rising up, 
setting up a "race condition" between the boundary and the walker. Lower panels D 
and E: 300 sample paths which start at the red point on the left and have their first 
passage through the boundary (white) on the red point in the right. White curve: 
average trajectory (analytic). Sample paths are colored by the probability density of 
the point they go through. In (D), hitting a left flank of a minimum is easy, and 
the average trajectory to do so does not significantly deviate from the deterministic 
trajectory until the very end, where the white curve can be seen to rise onto the 
minimum following a square root. In (E), hitting the right flank of a minimum is 
hard, and the average trajectory to do so strongly deviates from the deterministic 
trajectories of the system, missing the minimum by just enough not to collide with it, 
then rapidly rising to meet the first passage point, again, in a square-root profile. 



are generated by the walker running into a hard-to-avoid ob- 
stacle, as shown in Figure la. The intuition is that the walkers 
are moving left to right, rising onto a ceiling from which fea- 
tures are hanging, and as the walkers rise they collide with 
some feature. The problem is thus twice symmetry-broken: 
what matters are local minima of the boundary, not the max- 
ima, which are hard to get into; and the walkers only sponta- 
neously run onto the left flank of a local minimum. Therefore, 
a good first order approximation follows from observing that 
most of the first passages occur on the left flanks of local min- 
ima, and deeper local minima cast "shadows" on subsequent 
shallower minima. 

However, there is a finite probability that a walker may 
narrowly avoid a local minimum and pass just under it, only 
to rapidly rise afterwards and hit the right rising flank of the 
barrier, as shown in Figure 1C. This is, effectively, a race be- 
tween the boundary and the walker: if the walker can rise 
far faster than the boundary, then there is some probability 
of passage right of the minimum. But if the boundary rises 
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Fig. 2. Rasterplots and PSTH. A small segment of our dataset is displayed for 
clarity. A rasterplot and a plot of the PSTH are shown for each of three Holder 
exponents: 0.25 (rough), 0.5 (transition) and 0.75 (smoother, though still not differ- 
entiable). There's approximately the same number of spikes in all three groups. The 
rasterplots display the times at which the neuron fired (i.e. a first passage) stacked 
vertically (as a function of stimulus presentation number) to show repeatability. The 
PSTHs show a temporal histogram of said spikes. Please note the differences in ver- 
tical scale of the PSTHs: for Holder exponent H = 0.75 there are no bins with 
fewer counts than 10 events or more than 60, while for H = 0.25 most bins have 
counts while a few have over 1000 counts. 
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faster than a walker can catch up with, then the probability 
of passage right of the minimum can be exponentially small. 
Let us consider a local minimum of the barrier T(t) at time 
to of the form 

T(t > t ) « T(t ) + (t - t ) h 

and consider a walker that has just narrowly missed the min- 
imum by an amount e: W(to) = T(to) — e. The probability 
of the process to be at value W at time t > to is, to leading 
order, 



r(t - to) 

and thus the probability of arriving at the barrier at time t is 
approximately 



P h (t) = lim P(T(t),t) « exp(-(t ■ 

e— >0 



*o) 2H -Vr) 



When H < 1/2 this expression has an essential singularity, 
its value singular-exponent ially small for small times. In fact 
the probability and all of its derivatives are zero at to. For 
instance, consider a barrier whose flank to the right of the lo- 
cal minimum rises like v^At. As the fourth root in the barrier 
rises much more rapidly than the square root in the walker, 
the probability of hitting the barrier after the minimum looks 
like exp( — 1/y/At), a function that has an essential singularity 
at 0: the function as well as all of its derivatives approach 
as At -> + . 

The parameter H we described above, which is called the 
Holder exponent of the function, quantifies the ability of the 
barrier to, locally, rise faster or slower than a random walk. 
More formally, a function f(t) is said to be if- Holder con- 
tinuous if it satisfies \f(t) — f(s)\ < C\t — s\ h ; the roughness 




Fig. 3. (a) Probability density of firing as a function of time (horizontal) and 
Holder exponent (vertical), color coded in log scale. 51 values of the Holder exponent 
H between 0.25 and 0.75 are stacked vertically. The bin counts shown in the PSTHs 
of Fig 2 are color coded with a logarithmic code, (b) 3D rendering of a section of the 
data in (a): vertical axis and color scale is logarithmic in the rate, where it is evident 
that towards the back of the figure ( Holder exponent H = 0.25) the rate either 
diverges or goes to zero a.e. 



exponent H of the function is the largest possible value of H 
for which the function satisfies a Holder condition. Up to now 
we have considered a single local minimum, and even though 
the probability of crossing is singular-exponential small for 
H < 1/2, it is still nonzero. However, if the boundary is 
rugged, the local minima are dense. This density is not an 
issue for H > 1/2, that is inputs which are smoother than 
the internal noise; in this case the probability density of first 
passages is nowhere zero. But when H < 1/2 so the input 
is rougher or burstier than the internal noise, the probability 
density ceases to be a function: it is zero almost everywhere 
except for a set of zero measure where it diverges. 



Results 

We postpone to the Appendices the more formal proofs of reg- 
ularity of the first passage time probability distributions. We 
proceed now, instead, to discuss numerical simulations and 
their analysis. 

We carried out careful numerical integration of equation 
JTj|, for all Holder exponents H in the range (0.25 — 0.99) in 
increments of 0.01. In order for the results of the simulations 
at different Holder exponents to be directly comparable to 
one another, we generated the inputs I(t) by using the exact 
same overall coefficients in the basis functions of the Ornstein- 




Fig. 4. Density map of PSTH bin counts. The individual bin counts of the 
PSTHs as shown in Figs 2 and 3 are histogrammed here, and the value displayed 
as a logarithmic color map. All 7.5 billion spikes in our dataset were used for this 
plot. The bin counts are normalized by the average bin count (10 8 /2 22 ). For large 
Holder exponents, the probability of observing an actual count agrees with counting 
statistics given the average. As the Holder exponent becomes smaller, this distribu- 
tion becomes wider, until below 0.5 it becomes heavy-tailed. Notice the bottom row 
of the figure, representing the probability of observing a bin with zero counts. It is 
zero for all H > 0.5, becomes nonzero at H = 0.5, and for H < 0.5 it is the 
maximum of the distribution (i.e. the brightest red value). 



This transformation is referred as the Doob's transform. 
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Fig. 5. The tail of the cumulative probability distribution of observing a given count in the PSTH becomes a stretched exponential at Holder exponent H = 0.5. Top, 
the tails of the cumulative probability distribution, plotted as 1 — F(x) vs. x, for Holder exponents 0.4, 0.45, 0.5, 0.55 and 0.6 (right to left). The probability distribution 
is minus the derivative of these curves. Superposed on the data (black) a fit to the last 10 5 data points in the cumulative, i.e., the higher 2% percentile (red), in the form 
— log(l — F(x)) ~ ax + by/x + c. Right, the coefficients a, b and c for the aforementioned fit, plotted as a function of the Holder exponent H. Notice that the 
linear component a is (numerically) zero for H < 0.5, exposing the yfx term as the next higher order. For H > 0.5 the positive linear term guarantees convergence of all 
moments of the distribution. 



Uhlenbeck process described in [60], but scaled differently ac- 
cording to the Holder exponent laws (see Appendix). For 
each one of the 75 Holder exponents between 0.25 and 0.99, 
62000 repetitions of the stimulus were performed, accumulat- 
ing 100.000.000 first passages per Holder exponent. We com- 
puted the first passages using the fast algorithm described in 
[551 160] , which carries out exact integration in intervals which 
are recursively subdivided when the probability that the pro- 
cess attains the first passage exceeds a threshold, in our case 
10" 20 . The first passages were computed to an accuracy of 
2 -26 = 1/67108864, and the allowable probability that a com- 
puted passage is not in fact the first one is Pf a u = 10 -15 , so 
as to have an overall probability of 10 -5 that any one of our 
7.5 billion numbers is not in fact a true first passage. The 
values of the first passages were histogrammed in 2 22 bins; 
this histogram, which we call our PSTH (peristimulus time 
histogram) in analogy to the term in use in neural coding, 
represents the instantaneous probability distribution of first 
passage integrated over the bins, or, equivalent ly, the finite 
differences over a grid of the cumulative probability distribu- 
tion function for firing. 

The transition from smooth probability distribution to a 
singular measure is illustrated in Figures 2 and 3, where, as 
the Holder exponent is lowered, the concentration of the first 
passage probability on a small set is evident. Histogramming 
the individual bins of the PSTH we get the probability distri- 
bution to observe a given instantaneous rate of firing, shown in 
Figure 4. For large Holder exponents the rate does not deviate 
far from its mean. However, as the Holder exponent becomes 
1/2, both the probability of observing a zero rate, as well 
as the probability of seeing a rate far larger than the mean, 
become substantial. For H < 1/2 it becomes very probable 
to observe either zeros or large values of the instantaneous 
rate. This statement can be made precise by observing the 
tails of the probability distribution, and this is best accom- 
plished, given our numerical setup, by looking at the tails of 
the cumulative probability distribution, namely 




and then analyzing 1 — F{x) vs x for large x, which is carried 
out in Figure 5. Figure 5a shows that the tails of the distri- 



bution, when x ^> 1, decay exponentially for H > 1/2 but 
behave like stretched exponentials when H < 1/2: 

e~ ax , h>l/2, [2] 

l-F(x)^ e~ b ^ , h<l/2. [3] 

This observation is quantified in Fig 5b, where log(l — F) is 
fitted with a quadratic polynomial in namely 

— log(l — F(x)) ^ax + by/x + c 

For if < 1/2 the quadratic coefficient in the fit, which gives 
the convergent linear term, vanishes, uncovering the stretched 
exponential behavior. This quantitatively proves our assertion 
of a phase transition at H — 1/2. 

Discussion 

In abstract, mathematical terms, we have shown that the 
probability of observing a first-passage of a Gauss-Markov 
process through a rough boundary of Holder exponent H suf- 
fers a phase transition at H — 1/2. The integral of the prob- 
ability on equispaced grids becomes a stretched exponential, 
showing the underlying instantaneous probability has ceased 
to be a function: it is concentrated on a Cantor-like set within 
which it is infinite, and it is zero outside this set. Gauss- 
Markov processes, such as the Ornstein-Uhlenbeck process, 
can be mapped to the canonical Wiener process through a de- 
terministic joint scaling and time-change operation that pre- 
serves Holder continuity Furthermore, being the solution 
to a linear Langevin equation, the first-passage problem for 
drifted Gauss-Markov processes can always be formulated in 
terms of a fluctuating effective barrier that integrates the drift 
contribution. Therefore, our analysis directly applies to this 
situation. As non-linear diffusions with bounded drift behave 
like Brownian motion at vanishingly small scales, we envision 
that our result is valid for this more general class of stochas- 
tic processes with Holder continuous barrier. However, in this 
case, the barrier under consideration does not summarize the 
drift contribution of the diffusion. 

In terms of the original motivating problem, the encoding 
of an input into the timing of action potentials by a model 
neuron, this means that within our (theoretical and rather 
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aseptic) model, there is an abrupt transition in character of 
the PSTH, the instantaneous firing rate constructed from his- 
togramming repetitions of the same stimulus. The transition 
happens when the input has the roughness of white noise, con- 
ceptually the case in which the neuron is receiving a barrage 
of statistically independent excitatory and inhibitory inputs, 
each with a random, Poisson character. For inputs which 
are smoother than this, the PSTH is a well-behaved func- 
tion whose finite resolution approximations converge nicely 
and properly to finite values. However, when the input is 
rougher than uncorrelated excitation and inhibition, for ex- 
ample when excitatory and inhibitory activities are clustered 
positively with themselves and negatively with one another, 
then the PSTH is concentrated on a singularly small set, which 
means that the PSTH consists of a large number of sharply- 
defined peaks of many different amplitudes, but each one of 
them having precisely zero width. The width of the peaks 
is zero regardless of the amplitude of the internal noise; in- 
creasing internal noise only leads to power from the tall peaks 
being transferred to lower peaks, but all peaks stay zero width. 
Since the set of peaks is dense, refining the bins over which 
the PSTH is histogrammed leads to divergencies. 

Concentration of the input into rougher temporal patterns 
would evidently be a function of the circuit organization. For 
example, in primary auditory cortex, the temporal precision 
observed in neuronal responses [61 appears to originate in 
the concentration of excitatory input into sharp "bump" -like 
features [62] . 

It currently remains to be seen whether our mechanism 
will resist the multiple layers of real-world detail separating 
the abstract equation from real neurons in a living brain. 
Obviously, the infinite-sharpness of our mathematical result 
shall not withstand many relevant perturbations, which will 
broaden our zero- width peaks into finite thickness. That this 
will happen is indeed sure, but not necessarily relevant, be- 
cause a defining characteristic of phase transitions is that their 
presence affects the parameter space around them even under 
strong perturbations: that is why studying phase transitions 
in abstract, schematic models has been fruitful. Thus the 
real question remaining is whether our mechanism can retain 
enough temporal accuracy to be relevant to understand the 
organization of high-temporal-accuracy systems such as the 
auditory pathways, and whether our description of the rough- 
ness of the input as the primary determinant of coding modal- 
ity, temporal code or rate code, may illuminate and inform 
further studies. 



Appendix: Proofs 

Consider the stochastic leaky integrate-and-fire model for 
a spike triggering membrane threshold / and a post-spiking 
reset value r < I. Suppose a spike is emitted at time U > 0. 
With initial condition X t + = r, the inhomogeneous linear 
stochastic differential system 



dX t = -aX t dt + a dW t + dC(t) , t > U . 



[4] 



describes the ensuing sub-threshold noisy dynamic of the po- 
tential when driven by the input current dC(t). Here, dC(t) 
shall be considered as the infinitesimal increment of a time- 
varying load function C(t) that is if -continuous for a given 
Holder exponent H > 0, i.e. for every T > 0, there exists a 
constant ct > such that for all < t, s < T 



lim sup 

5^0+ | t _ a |<< 



\C(t)-C(s)\ 



< ct • 



Notice that, at the cost of rescaling X and I by a, we can 
restrain ourselves to the study of the case a = 1. 



Effective Barrier Formulation 

The nonlinearity of the leaky integrate-and-fire model lies en- 
tirely in the spike generation and subsequent reset, so that 
we can separately integrate input and noise between spikes. 
Thus, our first-passage problem for constant threshold / and 
varying forcing dC becomes a first-passage problem without 
drivingforces to a fluctuating effective barrier. Precisely, we 
solve Ejl writing X = U l -\-l l , where we separate the stochastic 
part JJ % (the Ornstein-Uhlenbeck process obtained for dC = 0) 
and the deterministic part P arising from the integration of 
the input dC(t): 



m 
i\t) 



re -a(t- ti ) + j\-«(t-s) 

^V**-^ dC(s). 



Determining the next spiking time U+i can be cast in terms 
of a first-passage problem for the process U l with the effective 
barrier t ^ L l {t) = I - l l (t): 



t i+ i = M{t > ti\U\ > V{t)}. 



[5] 



Therefore, a train of spikes to < t± < . . . < t n is determined 
by solving consecutively the first-passage problems [jS]] . Note 
that, due to the reset rule, the effective barriers do not agree 
at spiking times I? -1 (£~) / L x (tf) — I. However, for all 
i > 0, we have for t > U: 

{Ui<l-l\t)} = { u% t- £ e~ a{t ~ s) dC{s) </-f- 1 (t)} 



Making the left-hand term JJ't of the second inequality ex- 
plicit, we have 



v' t % 



-Oi{t-ti) 



rti 
Jti-i 



dC(s 



dW s 



and we recognize JJ't as the solution of Q for dC — 0, with 
the new initial condition: 



U t .+ =r- I e 

J ti-! 



} dC(s) 



L (t)-(/-r). 



As a result, the train of spikes to < t\ < . . . < t n is determined 
by the sequence of first-passage problem: 



U+i =inf{t > ti| if >L°(t)}, 



[6] 



where JJ' % is the standard Ornstein-Uhlenbeck process with 
initial condition U' t .+ — L°(t) — (I — r). In other words, by 
altering the reset rule, the linearity of the stochastic dynamics 
allows us to recast the successive first-passage problems pSI] 
in terms of a sequence of first-passage problems for one single 
continuous barrier L = L° f|6l] . 



First- Passage Markov Chain 

In a typical experiment, the spiking history of a neuron is 
recorded in response to repeated presentations of the same 
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stimulus. We idealize this situation by studying the distribu- 
tion of spiking events when an input cyclically forces a leaky- 
integrate- and-fire neuron. To avoid discontinuity effects, we 
choose a barrier satisfying L(T) = L(0) for some T > and 
then extend the definition of L on the whole time-line by peri- 
odization L(t) — L(t mod T). Then, the sequence of random 
times T n = (r n mod T), where r n denotes successive first- 
passage times to L, defines a discrete-time Markov chain T 
over the finite time period [0,T), seen as an oriented circkj^] 
To make it more formal, assume we can choose a load function 
satisfying for some T > 

f T e~ a(T - s) dC(s) = 0, [7] 
Jo 

which amounts to having a periodic effective barrier by set- 
ting L(t) = L(t mod T). For any time s in [0, T), consider the 
first passage time r s for an Ornstein-Uhlenbeck process start- 
ing at Us = L(i) — (I — r) and the barrier L. Because L(t) is a 
continuous function, it is known that the random variable r s 
admits a continuous non-decreasing cumulative distribution 
function F s : [s, oo) — > [0, 1] 70 . We then define the measure 
k s on the Borel sets of [s, oo) by setting for every open set 
O a ,b = (a, b) C [s, oo), s < a < b: 

k s (O a , b ) = F a (b) - F a (a) . 

Moving forward, we identify [0, T) with the circle C = R/TZ, 
which is compact for the Euclidean distance and for which the 
open arc circles 0( 0j b), are oriented counter-clockwise from a 
to 6, and generate the collection of Borel sets B(C). Equipped 
with the quotient map tt : R + ~ C, we define on the compact 
measurable state space (C, 8(C)) the measure kernels kj by 
setting for all open 0( a ,b) 

fc^(0( a , b )) = ks(7v~ 1 (0 (aib) )) . 

The collection of measures kj form a transition kernel on the 
compact state space C. Given an initial probability measure 
/io on C, they define a continuous state, discrete time Markov 
chain [671 [731 [75] T = (T,V) on (Sl,M) = (C, 23(C)) N , whose 
probability V satisfies: 

Vn e N, V(dr ni ... 1 dro) = 

k Tn _ 1 (dr n ) . . . k TO (dn)fio(dro) . 

In particular, for all u, v in C, v i— >• kj '(0( UjV )) is continuous 
in v with kJ(C) = 1. 

We shall see k s as the cumulative distribution of r n when the 
underlying process U n starts at U™ = L(s) — (I — r), i.e. the 
distribution of a spiking event knowing that the previous spike 
occurs at t. As such, the kernels kj need not admit a density 
k satisfying kj \dt) — k(s, t) dt, similarly to the "Devil's stair- 
case" resulting from the integration of the uniform measure 
over the triadic Cantor set [72] . 

Ergodicity of the Markov Chain 

We are interested in using this Markov framework to eluci- 
date the distribution of spiking events when a neuron is driven 
cyclically by an input defined [m] . To ensure that the instan- 
taneous firing rate and the probability of spiking coincide, we 
show that the Markov Chain (T, V) is ergodic, a notion we 
define in the following. 

An distribution \i is invariant by (T, V) if it satisfies 

f T T 

ji(dt) = k s (dt)fi(ds) , 



so that if Tn is distributed according to /x, so is Tn+i- When 
there exists a unique such measure /i, for any initial distribu- 
tion /io and any measurable set B on the circle C 

N — l , 

JyjEw^B)- !*(*) = { J If x 4b . 

and the Markov chain is said to be ergodic. Simply stated, 
the mean sojourn-time of the Markov chain in B tends toward 
the measure of B under ji. 

We can show that the Markov chain (T, V) is indeed ergodic 
for if -continuous functions with H > 0. Since the state space 
C of (T,V) is compact, it is enough to show that it has the 
strong Feller property [68] to prove the existence of invariant 
measures, i.e. 

V£e£(C), s n ^s£(C, k Sn {B) -> ks{B) . 

To establish the unicity of the invariant measure /i, it is 
enough to show that the Markov chain (T, V) has the irre- 
ducible property [68] : 

VBGi3(C), VsgC, k 8 (B)>0. 

We deduce the two properties above from consideration about 
the first-passage time problem in Supplementary Materials. 
The Feller property specifies that, if two identical leaky 
integrate-and-fire neurons spike respectively at times s and 
t, then, when s asymptotically approaches t, the probabil- 
ity that the first neuron later spikes in a given time interval 
becomes the same as for the other neuron. In other words, 
close initial conditions entail similar probability laws for the 
occurrence of the next spiking events (in the sense of the Kol- 
mogorov test). 

The irreducible property, which states that if one spiking time 
is achievable for a given starting condition (previous reset 
time), it is attainable for any starting time, similarly stems 
from these two intuitive observations. If one trajectory start- 
ing at t has a non-zero probability to hit this barrier in a given 
time region, we can easily convince ourselves that another tra- 
jectory starting at any s has a non-zero probability to be close 
to the reset value in t, and from there, unfold as a trajectory 
that has been reset in t. 

Intuitively, these properties holds for our first-passage Markov 
chain for two reasons. First, the continuity of the barrier 
which ensures the continuity of the cumulative distributions of 
the transition kernels. Second, the non-zero reset rules which 
constrain the membrane potential to be reset away from the 
barrier, thus avoiding pathological situations such as immedi- 
ate absorption. 

Numerical Simulation of the Markov Chain 

If the first-passage Markov chain (T, V) is ergodic, due to the 
possible irregularity of the barrier, numerical simulation of 
its invariant measure demands that we resort to an approxi- 
mation scheme. To justify this approach, we adapt a general 
result from [69 , clarifying in which sense a sequence of Markov 
chains (T N ,V N ) converge toward a limit chain (T,V) when 
N tends to infinity. 

Theorem 2 (adapted from [69 ): Let {X N , Q N ) be a sequence 
of strongly Feller Markov chains defined on a compact state 
space S. If, for any s in E>, the kernel probability measures 
of X N converge in law toward a limit probability measure q s , 
then, any limit in law of a sequence v n of invariant measures 
of (X N , Q N ), is an invariant measure of the Markov chain 



2 The passage of time orients the circle and we identify the future time T with the past time 
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(X, 2) corresponding to the limit kernel q. 

In particular, if all (X N Q N ) and (X ', Q) are ergodic, the 
sequence v n is uniquely denned and so is its limit distribution 
is, which is the stationary measure of (X, Q). 
For our purpose, an efficient approximation strategy of \i con- 
sists in exhibiting a sequence of ergodic strongly Feller Markov 



chains (T N ,V N ) 



whose kernels k 1 ^ converge to kj in law. 



This is accomplished by considering a sequence of first-passage 
Markov chains (T N , V N ) defined for the piecewise continuous 
periodic barriers Ln that interpolates L on the dyadic points 
D N {k2~ N T\0< k < 2 N }: 



L N : t e C \-> E \u t 



L(k2 



~ N T) ,0<k<2 N 



where E denotes the expectation with respect to the law of U 
(see [77]). Such Markov chains are ergodic by the same argu- 
ment as for (T,V). Moreover, since we restrain ourselves to 
barriers L that are H- continuous, the sequence Ln converges 
uniformly toward L (see Supplementary Materials), which in 
turn, implies the convergence in law (and in distribution) of 
kf toward kj . This demonstrates the cogency of approximat- 
ing L by L N . 



Frozen Noise as Injected Current 

In addition to providing a valid numerical method, the pre- 
vious approach provides an easy description of the input dC 
that gives rise to L. The central results is adapted from [78] : 

Theorem: There exists a Schauder basis of continuous func- 
tions ip n k compactly supported on S n h 

= [k2~ n+1 T, (k + 

l)2~ n+1 T] such that, for all N > 0, 

E \U t | U k2 -N T , < k < 2"] = ^n,k(t)^n,k 

0<n<N 0<k<2 n ~ 1 

where the ^ n ,k are the independent standard Gaussian vari- 
ables 

S n ,fc — / 4>n,k{t) dW t , (j)n,k — 1p'n,k ~ Ct1pn,k • 
J 

and the thus-defined functions n ,fc form an orthonormal sys- 
tem of L 2 [0, T]. 

Equipped with this result, it is easy to see that writing 
the input dC as a "Gaussian white noise" 

dC t = ^2 ^2 <t>n,k(t) -En,k , S n , fc i.i.d ~ JV(0, 1) , 

0<n 0<k<2 n - 1 

the statistics of the resulting random barrier 

L t = l- T e- a{t - s) dC{s) = 1-^2 J2 ' E ^ k ' 

^° 0<n <fe<2^-l 

is the same as for an Ornstein-Uhlenbeck process centered 
around zero and translated upward by /. Moreover, set- 
ting £o,o = 0, we naturally enforce the periodic condition 
L{t) = L(T) = I. 

However, we aim at studying the distribution of spiking events 
of a neuron cyclically driven by a deterministic input. Ac- 
cordingly, suppose now dC(t) — dCt(uj) is a realization of our 
"Gaussian white noise", i.e. a frozen noise. Then, L(t) = 
L t (uj) is the sample path of an Ornstein-Uhlenbeck bridge 
translated upward, which is almost surely if -continuous of ex- 
ponent 1/2. For this reason, we denote such an input dC 1 ^ 2 , 



Family of Holder Continuous Barriers 

From there, let us consider the set of coefficients ^ n ,k for 
which the continuous barriers of the form 

L N (t)=l~ ^2 ^2 VVfc(t) • £n,k , 

0<n<N 0<k<2 n ~ 1 

converge uniformly on C. It can be shown [78] that con- 
tains the set 

H£ = gR n |3 5 < 1,3 N > 0,Vn > iV,max|e n , fc | < 2 n<5/2 } 

k 

From this, we deduce that given L 1 / 2 , for any real H such 
that < H < 1, the barrier L H 



L H (t)=l~J2 J2 Mt)<n,k, Zn,k=2 
0<n 0<k<2 n ~ 1 



,n(ff-l/2) 



£,n,k 5 



is well-defined as a continuous function of C. Keeping this in 
mind, we have at our disposal a well-known result [65] relating 
the local Holder exponent of a function to the asymptotic be- 
havior of the coefficients of its decomposition in the Schauder 
basis. Adjusting to our situation, it directly entails that for 
all if, < H < 1, the barriers L H are almost-surely H- 
continuous. Therefore, we can continuously (in the L°°-norm) 
control the asymptotic Holder continuity of the effective bar- 
rier driving the activity of a leaky integrate-and-fire neuron 
by smoothly changing the coefficient £ n?fe used to construct 
piecewise approximations L%. 

In order to emphasize the effect of the varying Holder reg- 
ularity, we adopt a slightly modified version of our barriers 
L H , by weighting them with a continuous function H H> c{H) 
under the from L' H = c(H) (L H - L H (0)) + L H (0). The func- 
tion c is chosen so that the newly formed barriers cause the 
neuron to fire with an overall mean firing rate (as opposed to 
the instantaneous mean firing rate which is time-dependent) 
remains constant when changing H. Formally, this constraint 
is equivalent to holding a constant mean inter-spike time 



/T / poo 
(J (t-s) K f(dt)) fi H {da) 

while varying L^ 



Integral Equation for the First-Passage Time 

We establish the existence of a density function for the first- 
passage time of a Wiener process hitting a if-continuous bar- 
rier with H > 1/2. This property is formally referred to as the 
absolute continuity of the first-passage time distribution with 
respect to the Lebesgue measure on the real half-line. With- 
out loss of generality, we adopt the point of view of a killed 
Wiener process absorbed on a fluctuating boundary, which al- 
lows us to use the powerful machinery of the heat equation. 
The presented result stems from the ground-breaking work of 
Gevrey [66] about parabolic differential equations, later actu- 
alized in a modern form by Rozier [63] . 

Integral equations for the cumulative distribution of the first- 
passage time of a Wiener process naturally arise from proba- 
bilistic arguments. Consider the event {W t > x} for an con- 
tinuous barrier L satisfying x > L(t). Then, the first-passage 



the associated barrier L 1 ^ 2 and the coefficients 



1/2 



Notice that for the sake of well-posedness, the kernels that intervene in the formulation of the 
mean inter-spike time are computed for a periodic barrier L H but defined on [0, oo) instead of 
being wrapped on [0, T). 
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time r with L occurs certainly before t and we can condition 
this event with respect to r, which yields 



P(Wt > x) = E [¥(W t >x\t) 



[ F(W t 
Jo 



> x I r — s)q(ds) 



[8] 

where q denotes the first-passage time probability measure. 
Using the strong Markov property, on {r = s}, we can dis- 
regard the past-trajectory of W and equate the probabilities 
¥(W t > x\r = s) and F(W t - s > x - L(s)). Differentiating 
equation J8| with respect to x, we end up with 



x 

~7t 



Jo 



L(s) 



q(s) ds , 



where k denotes the Heat kernel. 

It is important to observe that as long as L is if-continuous 
with H > 1/2, we have 



lim 



Lit) ~ L(t) 

Vt^ 



= 0. 



Since k is a smooth function, we can make the arbitrary value 
x tend toward the barrier L(t) by superior value and, through 
the dominated convergence theorem, we get the following in- 
tegral Volterra equation [?, ?]: 



L{t) - L(s) 



q(s) ds . 



[91 



This integral equation, which dates back original work from 
Siegert [74], stems from the fact that k(s,x; W s ,t) indexed by 
s is a martingale [76], which offers a convenient way to gen- 
eralize this equation to general time-inhomogeneous diffusion 
processes. 



Absolute Continuity of the First-Passage Time 

The integral equation is of the Volterra type, which comes 
in two flavor: equations of the first kind and of the second 
kind 71 . To ensure the existence and unicity of a solution 
to the equations of the second kind, we have the following 
powerful result: 



Theorem (adapted from [631IT9] ) 
tion of the second-kind 



The linear Volterra equa- 



9(t) =/(*)+ f K(t,s)f(s)ds, 
Jo 

where g is a piecewise continuous function has a unique piece- 
wise continuous solution / for all t > if K is bounded on 
< s < t and if there exists a monotone increasing function 
a with linit^o a(t) = 0, such that for all < s < t 



r. 



\K(t,r)\dT < a(t-s). 



Unfortunately, equation [p)] is a Volterra equation of the 
first-kind and as such cannot oe dealt with directly. However 
for barriers L that are i7-continuous, it can be recognized as a 
linear generalized Abel integral equation, that is an equation 
of the type 



git) 



-r 

J s 



(t-r) h 



dr 



where / is the unknown, g is a continuous function, and K is 
a continuous kernel for s < t and < h < 1. 
Abel integral equations are frequently encountered in physics 
and there are methods to prove the existence and unicity of a 
solution by transforming the original equation into a Volterra 
equation of the second-kind. In our case, it proceeds through 
the use of the Abel integral transform, which is designed to 
solve the canonical Abel equation 



g(t) 



f 

J s 



fir) 



The unique solution is given as 



dr. 



f(t) = A[g]{t) = 



7i dt 



f 

J s 



9(r) 



dr 



where A is the Abel inverse operator. The application of A 
to equation [JoJ reduces the problem to a Volterra equation of 
the second-kind: 

Proposition (adapted from 63 ): If L is if -continuous with 
H > 1/2, through the application of the Abel operator, the 
Volterra equation of the first-kind Q is equivalent to the 
Volterra equation of the second-kind 



2tt A[g] (t) = q(t) + - f K(t, r)q(r) dr , 

* Js 

with the kernel K being defined as 

(l(ct)-L(t)) 2 



K(t,T) 



d_ 
dt 



y/(t-a)(a-r) 



dr [ 



and g denotes the continuous function g(t) = fc(s, x\ t, L(t)) . 

A careful study shows that the kernel K satisfies the con- 
ditions of Theorem [63]. Thus the integral equation [10] 
obtained through the Abel transform admits a unique contin- 
uous solution, which is the density of the first-passage time to 
the barrier L. 
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